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We show that classical molecular density functional theory (MDFT), here in the homogeneous reference fluid 
approximation in which the functional is inferred from the properties of the bulk solvent, is a powerful new 
tool to study, at a fully molecular level, the solvation of complex surfaces and interfaces by polar solvents. 
This implicit solvent method allows for the determination of structural, oricntational and energetic solvation 
properties that are on a par with all-atom molecular simulations performed for the same system, while reducing 
the computer time by two orders of magnitude. This is illustrated by the study of an atomistically-resolved 
clay surface composed of over a thousand atoms wetted by a molecular dipolar solvent. The high numerical 
efficiency of the method is exploited to carry a systematic analysis of the electrostatic and non-electrostatic 
components of the surface-solvent interaction within the popular CLAYFF force field. Solvent energetics and 
structure are found to depend weakly upon the atomic charges distribution of the clay surface, even for a 
rather polar solvent. We conclude on the consequences of such findings for force-field development. 

PACS numbers: 61.20.Gy, 61.20.Ja, 61.25.Em, 68.03.Hj, 68.08.-p, 68.08.De 



I. INTRODUCTION 

Liquids at solid-liquid interfaces behave differently 
from bulk liquids. Intcrfacial phenomena play key roles 
in various applications like, for instance, heterogeneous 
catalysis, electrochemistry, adsorption and transport in 
porous media. It is our goal to improve the description 
of such interfaces at the microscopic scale. 

Experimentally, specific techniques have been devel- 
oped, such as the vibrational sum frequency spec- 
troscopy^ or quasi-elastic neutron scattering techniques^, 
that now give access to information on the structural and 
dynamic properties of the liquid on few molecular layers 
on top of the surface. At these scales, theoretical model- 
ing can clarify experimental observations by linking them 
to atomistic phenomena. For instance, Rotenberg et al£ 
showed by molecular simulations that wetting properties 
of talc surfaces, that behave as hydrophilic or hydropho- 
bic depending on the relative humidity, arc a consequence 
of the competition between surface-solvent adhesion, fa- 
vorable entropy in the gas phase and molecular cohesion 
in the liquid phase. 

This kind of theoretical modeling is typically based 
on molecular dynamics (MD) or Monte Carlo simula- 
tions where each atom is considered individually. For 
systems that are intrinsically multiscale, e.g. porous me- 
dia, one has to deal with thousands of atoms. At this 
point, numerical heaviness becomes critical and hinders 
systematic studies. To overcome this difficulty, various 
mesoscale models have been proposed, that are less com- 
putationally demanding. One can for instance forget the 



a ' Electronic mail: maximilicn.lcvcsque@gmail.com 



molecular nature of the solvent to only account for a Po- 
larizablc Continuum Medium (PCM)^ - — , or simplify it to 
a hard sphere, with or without dipol© 7 - - — . Coarse-grained 
strategies have also been successfully applied to solid- 
liquid interface a 10 ' 11 . While these methods and others 
proved efficient on a panel of problems, the applicability 
of their underlying approximation to an unknown system 
remains to be carefully and systematically checked. 

In the second half of the last century, liquid state 
theorie s 12 ' 13 like the classical density functional theory 
(DFT)Mri6 have blossomed. Classical DFT, like other 
methods such as integral equations or classical field the- 
ories have shown to give thermodynamic and structural 
results comparable to all-atom simulations at attractive 
numerical cost, but they have for now been limited to 
highly symmetric systems, bulk fluids or simple model 
interfaces, e.g. hard and soft walls. A current challenge 
lies in the development of three-dimensional theories and 
implementations to describe molecular liquids, solutions, 
mixtures, in complex environments such as biomolecular 
media or atomically-resolved surfaces and interfaces. Re- 
cent developments in this direction use lattice fiel d 17 ' 18 , 
or Gaussian field^ theoretical approaches, or the 3D 
reference interaction site model (3D-RISM)22r— , an ap- 
pealing integral equation theory that has proven recently 
to be applicable to, e.g., structure prediction in com- 
plex biomolecular systems. Integral equations are, how- 
ever, restricted by the choice of a closure relation (typ- 
ically HyperNetted Chain (HNC), Percus-Yevicb^, or 
Kovalenko-Hirata 2 ^ approximation), and despite their 
great potential, they remain difficult to control and im- 
prove, and can prove difficult to converge. 

Recently, a molecular density functional theory 
(MDFT) approach to solvation has been introduced^ - — 
It relies on the definition of a free-energy functional de- 
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pending on the full six-dimensional position and orienta- 
tion solvent density. In the so-called homogeneous ref- 
erence fluid (HRF) approximation, the (unknown) func- 
tional can be inferred from the properties of the bulk 
solvent. Compared to reference molecular dynamics cal- 
culations, such approximation was shown to be accurate 
for polar, non-protic fluids 3 ^ - — , but to require correc- 
tions for wate r 31 i 33 ~ — . Until now, only simple three- 
dimensional solutes have been considered in this frame- 
work including ions and pola r 31 ' 32 or apolar— molecules. 
In refi22, the authors have shown that the MDFT scheme 
was applicable to the computation of electron-transfer 
properties such as reaction free energies and solvent re- 
organization free energies. In this paper, we show how 
the most recent developments of the molecular density 
functional theory, described here in the homogeneous ref- 
erence fluid approximation 2 ! (HRF-MDFT), can be im- 
plemented to study the solvation of a complex atomically- 
resolved clay surface, the pyrophyllite, that contains over 
a thousand atoms. In order to uncouple the difficulties 
intrinsic to the method to those due to the particular cor- 
relations in water, we restrict ourselves here to a generic 
model of dipolar solvent, the Stockmayer fluid - with 
eventually parameters adjusted to mimick a few proper- 
ties of water. 

In section HH we present the most recent developments 
in molecular density functional theory in the homoge- 
neous reference fluid approximation. We also describe 
our reference molecular dynamics simulations and the 
system under investigation. In section IIII1 we discuss 
the structural and orientational properties of the com- 
plex interface between the solid and the dipolar liquid. 
There, we compare HRF-MDFT with reference molecu- 
lar dynamics results. Finally, in section IIII CI we illus- 
trate the possibilities offered by the method in terms of 
numerical efficiency to analyze the electrostatic and non- 
electrostatic components of the state-of-the-art force field 
for clays, CLAYFF— . We conclude on the consequences 
of our analysis on force field development. 



II. METHODS AND MODEL SYSTEM 

A. Molecular Density Functional Theory of Solvation 

1. Theoretical Aspects 

Recent developments of a molecular density functional 
theory (MDFT) unlocked the study of the solvation of 
three-dimensional molecular solutes in arbitrary solvents 
using classical density functional theory. In MDFT, the 
solvent is composed of rigid molecules and described in 
terms of a position and orientation density, p (r, f2). The 
solvation free energy is defined as the difference between 
the grand potential of a system that includes exter- 
nal perturbation (a solute) and inhomogeneous solvent, 
and the grand potential &[po] of the homogeneous solvent 
at density po = no /8ir 2 , where no is the number density, 



and without external potential, 

F\p] = $\p]-*\po], 



(1) 



where p and po are the position and orientation densities 
of the inhomogeneous and homogeneous solvent. Follow- 
ing the theoretical framework introduced by Evans^— , 
the density functional F can be rewritten as the sum of 
an ideal contribution (Fid), an excess term (F exc ) and an 
external (solute) contribution (F ext ), 



F[p] = F id [p] + F ext [p] + F exc [p\. 



(2) 



The ideal and external contributions F,,j and F ex t can be 
formally and exactly expressed as 



F id [p]=k B T J drdfl[p (r,0)ln 

-p(r,fl) + po], 
F ext [p} = J drdfl V ext {r,n)p(r,n), 



(3) 
(4) 
(5) 



where k B is the Boltzmann constant and T the temper- 
ature. V ex t is the sum of the electrostatic and Lennard- 
Jones interactions between the solute and one solvent 
molecule located at r with orientation f2: 

V ext (r,Q)= ]T {qjV^rj) 

j £ solvent 



i £ solute 
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If R(Q) is the rotation matrix associated to il and Sj is 
the position of the solvent site j in the molecular frame, 
then Yj = r + R(f2)sj denotes its absolute position in 
space and ry = rj — r, its relative position with respect 
to the solute site i located at . e,j , ay are the Lennard- 
Jones parameters between solute site i and solvent site 
j. qj is the partial charge carried by site j and V q {jCj) is 
the electrostatic potential created by the solute at rj. 
The excess functional is unknown, but can be expressed 
formally as 

F exc [p] = -\k B T J j dridr 2 dn 1 dn 2 &p{ri,tti) 

xc(r 2 - ri,fi 1 ,0 2 ) Ap(r 2 ,Q, 2 ) 

+F B [Ap], (7) 

where Ap (r, fi) = p (r, fi) — p . The first term rep- 
resents the homogeneous reference fluid approximation 
where the excess free-energy density is written in terms 
of the angular-dependent direct correlation of the pure 
solvent. It amounts to a second order Taylor expansion 
of the excess free-energy functional around the homoge- 
neous solvent. It is equivalent to the HNC approximation 
in integral equation theories when the solute is taken as 
a solvent particle. The second term represents the un- 
known correction to that term (or bridge term) that can 
be expressed as of a systematic expansion of the solvent 
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correlations in terms of the three-body,. . . n-body terms 
direct correlation functions. We will consider below the 
case Fb = 0. This approximation was shown in Ref.— to 
be accurate for several polar aprotic solvents. Difficul- 
ties were found, however, for water, certainly the most 
interesting but also the most complex solvent, in part 
due to its high-order correlations, if not interactions. We 
have proposed several corrections for water entering in 
the definition of Fb ■ an empirical three-body correlation 
terrnjU j nS pi rec } by the water model by Molincro et aZ.— , 
and a bridge function extracted from a hard sphere func- 
tional^. The first one enforces the missing tetrahedral 
order, while the second one introduces the iV-body re- 
pulsion terms (N > 2) of a hard-sphere fluid. 
Since water adds some extra complexity in the func- 
tional that deserve to be tackled separately, it is not 
our purpose to study hydration here, but rather to show 
the advantages/disadvantages of MDFT and its practi- 
cal implementation for the molecular solvation of com- 
plex atomically-resolved surfaces. We will thus consider 
a generic dipolar fluid, the Stockmayer model, with pa- 
rameters adjusted to mimick a few properties of water 
(particle size, molecular dipole and dielectric constant), 
but no hydrogen bonds, and certainly not a "drinkable" 
water. 

For dipolar models, each orientation 57 can be described 
by the normalized orientation vector ft and the direct 
correlation function in Eq. [7] may be developed on a ba- 
sis of three rotational invariants 



c(ri2,r2i,fi 2 ) 



c S (r 12 )$ 100 

+c A (ri 2 )$ 110 (fii,n 2 ) 
+c D (ri 2 )* 112 (ni,n a ) ; 



where 



$ ioo = i 



$ 110 = fix n 2 , 



(8) 



(9) 
(10) 



$ 1 



3(ni-fia)(n 2 -ria)-ni-n a> (11) 



with ri2 = r 2 — i"i and f 12 the associated unit vector. 
If we also define the molecular density n (r) and the po- 
larization density P (r) 



n(r) 



p(r, ft) dfl, 



P(r) = J ft ■ p (r, ft) dfl, 



(12) 
(13) 



the excess free-energy density functional F exc [p] can 
be rewritten as a functional of n and P instead of 
p(r,Q)2L28 : 

F exc [p(r,ft)} = ~k B T J dridr 2 {cs(ri2) An(n) An(r 2 ) 

-CA(ri 2 )P(ri)-P(rx) (14) 

-ci ? (r 12 )[3(P(r 1 )-f 12 )(P(r 2 )-f 12 ) 

-P(r x )-P(n)]}. 



This significantly increases numerical efficiency, as will 
be discussed below. The same reduction is true for the 
ideal and external contributions if the solute-solvent elec- 
trostatic interaction is strictly restricted to charge-dipolc 
interaction s) 27 ! 28 . Even in that case, however, it remains 
advantageous, both for convergence reasons and for keep- 
ing the generality of the code, to stick to the expression 
of the ideal free-energy in terms of the angular distribu- 
tion p (r, ft), eq. |4j and to minimize the functional in the 
full position-angle space. This is now described. 



2. Numerical Aspects of HRF-MDFT 

Here, we give insight into numerical details associated 
with the variational minimization of the density func- 
tional F, i.e. the resolution of the Euler-Lagrange equa- 
tion 



SF[P] 
Sp 



ksT In 



Po 



-V ext ( r ,tt) + ^f^=0. 

Sp 

(15) 

The inhomogeneous density p is projected onto an or- 
thorhombic position grid of N x x N y x N z nodes with pe- 
riodic boundary conditions. To each node is associated 
an angular grid on which is discretized the orientation 
vector ft. The variational density p (r, ft) which mini- 
mizes F is optimized numerically by the limited memory 
Broydcn-Flctchcr-Goldfarb-Shanno (L-BFGS) method 
as implemented by Byrd, Lu, Nocedal and Zh u 38 i 39 . This 
quasi-Newton algorithm only requires the knowledge of 
the functional F and its first derivative with respect to 
the density at each node, which are known analytically. 
The Hessian, needed in Newton-derived algorithms, is 
approximated using gradients at previous self-consistent 
iterations. This makes a notable difference with other op- 
timizers: faster than Piccard iterations, without requir- 
ing the second variation of the functional with respect 
to the density as typically required by pseudo-Newton 
algorithms. Convolutions in F exc are calculated in re- 
ciprocal space by fast Fourier transforms (FFT) as im- 
plemented in the FFTW3 library^ - - — . Rewriting F as a 
functional of n and P reduces considerably the number of 
angular summations and of FFTs to be performed. As a 
summary, the high performance of our three-dimensional 
HRF-MDFT is due to (i) a quasi-Newton implementa- 
tion, (ii) the calculation of convolutions in reciprocal 
space, associated to the extreme performance of FFTW3, 
(Hi) the rewriting of the functional in a form optimal for 
our purpose. 

For the given solvent molecular model, the direct cor- 
relation function of the solvent is extracted from all-atom 
simulations from which are generated the pair distribu- 
tion functions. The corresponding correlation function 
can then be deduced by solving the Ornstcin-Zernikc 
equation. This can be done in Fourier space but raises 
numerical issues at large r, i.e. at small k, the conju- 
gate of r. The direct space method of Baxter combined 
with the variational method of Dixon and Hutchinson 
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FIG. 1. Projections cs (black line), ca (red dashed line) and 
cd (blue dotted line) of the direct correlation function of the 
Stockmayer fluid on the first three rotational invariants. 



was used accordingl y 27 ) 43 . This method enjoins that c 
vanishes beyond a radius, set to 8.7 A in the present sys- 
tem. The projection of the direct correlation function of 
the Stockmayer fluid on rotational invariants c$, ca and 
Cd as expressed in Eqs. [5] to QT] arc plotted as a function 
of k in figure [T] Note that the exact bridge function of 
the Stockmayer fluid from explicit Monte Carlo simula- 
tion data has been published recently by Puibasset and 
BelloniM. 

The external potential V exc (r, SI) is computed only 
once, at the beginning of the simulation: At each node, 
an isolated solvent molecule is considered with a given 
molecular orientation, for which is tabulated the total in- 
teraction energy between the solvent molecule sites and 
the solute sites, according to Eq. [51 To accelarate the 
computation, the value of the electrostatic potential at 
each solvent site is interpolated from its values on the 
grid. The grid electrosatic potential itself is obtained 
by first extrapolating the solute charge density on the 
grid and then solving the resulting Poisson equation by 
FFT's. The calculation of the Lennard- Jones external 
potential takes advantage of cut-off distances for each 
solute site. 

In this work, we used typically 4x4x4 nodes per 
A 3 , which, for the system considered below, corresponds 
to roughly 150 3 3D-grid points. The molecular orienta- 
tion SI was discretized over 18 angles and integrated over 
the whole sphere by Legendre quadratures. This implies 
4x4x4x18 = 1152 variables per A 3 to be optimized. All 
results presented here were carefully checked with respect 
to the number of nodes per A 3 and per orientation. The 
convergence of the results as a function of grid resolution 
can be quantified as follows: with respect to calculations 
with a very fine grid resolution of 5 3 grid points per A 3 , 
the relative difference in solvation free energy amounts to 
1.2 % with 2 3 grid points per A 3 and 0.3 % with 3 3 grid 
points per A 3 . It is less than 0.1 % with 4 3 grid points 
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FIG. 2. a) Variation of the solvation free energy functional 
F between two iterations normalized by the initial value Fo 
during the iteration process. The grid is composed of 4 x 4 x 4 
nodes per A 3 and 18 discrete orientations per node. A typical 
minimization is performed with a precision of 10 in the 
solvation energy in less than 15 iterations, b) CPU time in 
seconds per iteration for grid meshes N x x N y x N z = 2 3 , 3 3 , 
4 3 and 5 3 per A 3 and 18 discrete orientations per node. The 
calculation of the external potential for a given grid mesh is 
done only once at the beginning of the simulation and takes 
a CPU time comparable to 1 to 2 iterations. The solvated 
solute consists of 1280 point charges and 1152 Lennard- Jones 
sites, for which parameters are extracted from the popular 
CLAYFF force neld 3 ^ and given in table [U The supercell 
volume is approximately 68 nm 3 . 



per A 3 . Values given above are found to be adequate 
for all observables presented in this paper. In order to 
illustrate the high efficiency of the method, the iterative 
convergence of F with iteration steps is illustrated in fig- 
ure [2] for the grid described above. After 15 iterations, 
the convergence in solvation free energy is of the order of 
10 -5 . We also illustrate in figure[2]the CPU time per it- 
eration for 18 angles per node and 2 3 , 3 3 , 4 3 and 5 3 nodes 
per A 3 for our model system containing several hundred 
atoms. One observes a linear scaling in N x x N y x N z , 
which leads to convergence in less than 15 minutes for 
the 4x4x4 grid per A 3 described above on an ordinary 
laptop^, without parallclization. 



B. Molecular Dynamics 

Molecular dynamics simulations were generated for the 
same surface and solvent molecular models in order to 
compare with the MDFT results. All simulations were 
done in the canonical NVT ensemble, with a Nose-Hoover 
thermostat. After a phase of equilibration, all structural 
quantities were collected and averaged over 5 ns. The 
local molecular density n(r) of the Stockmayer fluid is 
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Top view 
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FIG. 3. Microscopic clay structure, side and top views (red, 
O; white, H; yellow, Si; cyan, Al). In the top view, only lay- 
ers above the Al atoms are shown to highlight the hexagonal 
symmetry of the first O-Si layer, the hydroxyl group (-OH) 
parallel to the sheet, and the Al position in 2/3 of the octa- 
hedral sites. 



Molecule 


Atom 


e (kJ/mol) 


a (A) 


3(e) 


Pyrophyllite 


Al 


5.56388e-6 


4.27120 


1.575 


Si 


7.7005e-6 


3.30203 


2.100 


O g 


0.650190 


3.16554 


-1.050 


Oh 


0.650190 


3.16554 


-0.950 


H G 


0.0 


0.0 


0.425 


Stockmayer 


central 


1.847 


3.024 


0.0 


side 1 


0.0 


0.0 


1.91 


side 2 


0.0 


0.0 


-1.91 



TABLE II. Force field used to model the pyrophyllite solute 
and the Stockmayer fluid in both molecular dynamics simu- 
lations and HRF-MDFT minimization. Pyrophyllite param- 
eters are extracted from the CLAYFF force field^S. 



type 





Si 


O 


H 


Al 


H 


O 


Si 


O 


*(A) 


0.00 


0.59 


2.18 


2.27 


3.27 


4.27 


4.36 


5.95 


6.54 



TABLE I. Coordinate z along the normal to the surface of 
each site of the pyrophyllite. The distance Az between two 
successive layers is also indicated. 



calculated as the time averaged density in elementary 
volumes of 0.1 3 A 3 . The densities perpendicularly to the 
clay layers can be calculated by averaging the above den- 
sity in the (x,y) planes. The spatial densities of the av- 
erage dipoles were calculated in the same way. Molecular 
dynamics simulations were performed with the DLPOLY 
packag o 46 ' 47 . 



C. System Description 

Pyrophyllite is neutral clay. It is monoclinic of space 
group 2/m, perfectly cleaved along orientation {001}. 
Two views of a pyrophyllite sheet are provided in fig- 
ure [3] It consists of a stack of 6 atomic layers. The top 
layer consists of oxygen (O) and silicon (Si) atoms with a 
lateral hexagonal symmetry. The central layer consists of 
aluminum (Al) atoms in 2/3 of the octahedral sites. Be- 
tween these two layers is a oxygen-hydrogen (O-H) layer, 
with O atoms at the center of the hexagons formed by the 
top Al-0 sites. The O-H axis is oriented in the direction 
of the empty octahedral site. The coordinate z of each 
site along the normal to the clay surface is given in tableU 
The simulation box contains two half clay layers of lat- 
eral dimensions L x x L y = 41.44 x 35.88 A 2 , which corre- 
sponds to 32 clay unit cells of formula ALJSisC^oKOH).!. 
The distance between the surfaces is L z = 45.57 A, cho- 
sen to recover the bulk density of the Stockmayer fluid 
at the center of the pore, i.e. n = 0.0289 molecule per 
A 3 . The simulation supercell is thus w 68 nm in peri- 
odic boundary conditions, and contains 1280 clay atoms 
(640 per half clay layer). This pore contains 1600 solvent 
molecules in the reference all-atom simulations. 

We used the CLAYFF force field for both molecular 



dynamics simulations and MDFT minimization. It is a 
general purpose force field for simulations involving mul- 
ticomponent mineral systems and their interfaces with 
solutions. Related information is tabulated in table [TTJ 
The Stockmayer fluid is described by a 3-sites molecule. 
The neutral central site interacts via Lennard- Jones po- 
tentials only. The external sites have opposite charges 
of ±1.91 e, both distant of 0.1 A from the central site. 
It results in a solvent molecular dipole of 1.84 D (inci- 
dentally, approximately that of a water molecule in the 
gas phase), and, with the chosen Lennard- Jones param- 
eters that match the reduced parameter set of Pollock 
and Alder—, in a dielectric constant of roughly 80, com- 
parable to that of bulk water at room conditions^ 7 -. The 
parameters for the Stockmayer fluid force field can also be 
found in table [Hj Once again, it is not our purpose here 
to study the more complex solvation by water, which re- 
quires extra terms in the functional, but to demonstrate 
the possibilities of MDFT in the HRF approximation for 
a generic polar solvent whose functional is of a good qual- 
ity. 



III. RESULTS 

We first compare the predictions of HRF-MDFT for 
the solvent density and orientation to reference all-atom 
simulations, which are two orders of magnitude slower to 
acquire. We then analyze the role of electrostatic inter- 
actions on these quantities. 



A. Density Profiles 

The main structural observable describing solvated in- 
terfaces is the number density profile n(z), defined as the 
average of the number density on plane z, 

n{z) = — !- / / ^-dxdy. (16) 
L> x L y J J n 
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z(A) 

FIG. 4. Normalized number density profile of the solvent 
between two pyrophyllite layers as calculated by molecular 
dynamics (open circles) and HRF-MDFT (black line). 



The density profiles extracted from explicit molecular dy- 
namics and from HRF-MDFT are given in figurelH Good 
overall agreement is found between the two methods. A 
shoulder followed by two main peaks are found in MD. In 
HRF-MDFT, the shoulder looks more like a weak peak. 
The two strongest peaks are followed by long-ranged os- 
cillations. The first feature (shoulder or weak peak) is 
found at z = 8.3 A, 1.76 A away of the surface layer 
(see Table U to identify layers coordinates) . Its inten- 
sity is limited to w 0.7 in MD and « 1 in HRF-MDFT. 
The largest and main peak is at z — 9.5 A, 2.96 A af- 
ter the top surface layer. An intermediately intense and 
broad peak is also found at z = 12.3 A, 5.76 A away of 
the surface. These three structures will later be called 
"prepeak" , "main peak" and "secondary peak" . Further 
weak oscillations are found up to 15 A away from the sur- 
face. At the center of the supercell, the density is flat at 
n = no , namely at the bulk density, which means for sol- 
vent molecules to be in bulk, homogeneous, conditions. 

The only noticeable difference lies in the weak prepeak 
seen in HRF-MDFT observed as a shoulder in molecular 
dynamics. The HRF-MDFT is known to slightly overes- 
timate the height of the first peak in polarized systems^!. 
The localized nature of the prepeak can be seen in figure[5] 
where is presented the density map in the plane defined 
by the maximum of the prepeak, as calculated by MD 
and HRF-MDFT. The same position and overall shape 
is found with the two methods. It is slightly broader 
in-plane in HRF-MDFT for approximately the same in- 
tensity, which induces the higher value once averaged in 
the plane. The prepeak is localized at the center of the 
hexagons formed by surface Si and O atoms. One may 
note the high maximum value of the normalized number 
density n/n there (up to s=s 30). The integral of the 
density in this peak is the total number of particles in 




FIG. 5. Maps of the solvent normalized number density 
n(r)/no in planes of the prepeak (left), main peak (center) 
and secondary peak (right), as calculated by molecular dy- 
namics (top) and HRF-MDFT (bottom). 

it. While a strict deconvolution of the three-dimensional 
prepeak is arbitrary, no consistent definition results in 
more than one solvent molecule inside each prepeak. 

In figure[SJ we also plot the number density in the plane 
of the main and secondary peaks. Again, the shapes are 
very similar in MD and HRF-MDFT. The main peak is 
localized on top of Si atoms. On the contrary, a depletion 
is found on top of O atoms of the surface layer. The broad 
secondary peak can be found again on top of the center 
of hexagons, i.e. on top of O atoms. 

We now have a clear three-dimensional view of the sol- 
vent structure on the pyrophyllite surface: rare solvent 
molecules are adsorbed very close to the surface, at the 
center of hexagons formed by Si and O atoms of the clay 
surface layer. On top of these molecules is stacked a 
strongly structured hexagonal layer of solvent molecules, 
above Si atoms. An additional, more diffuse layer is 
found once again on top of the center of the hexagons. 
The relatively small distance between the first two lay- 
ers demonstrates significant interactions, and thus cohe- 
sion, between layers. This conclusion is in agreement 
with Rotcnbcrg et al. who showed recently how the com- 
petition between adhesion and cohesion determines the 
hydrophobicity of these surfaces^. 

As a partial conclusion about structural properties, HRF- 
MDFT results are in quantitative agreement with the ref- 
erence molecular dynamics simulations, with a speed-up 
of two orders of magnitude. 

B. Orientational Properties 

We now turn to the orientational properties of the 
molecules solvating the pyrophyllite surface. Even if it 
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is routinely done, it is worth emphasizing that orienta- 
tional properties are slower to sample in molecular dy- 
namics or Monte Carlo simulations as each elementary 
volume described in section Hi Bl has to be sampled for a 
number of angles. This problem is recurrent in numer- 
ous simulation techniques where orientational degrees of 
freedom (from electronic spins to molecular orientations) 
induces new and rich behaviors^—. For its part, HRF- 
MDFT gives a direct access to this quantity through the 
full orientational density p(r, fl), and also the polariza- 
tion density P(r) defined in Eq. [T31 as one of the two 
natural observables of the theory. 

The polarization density is found to be aligned in the 
z direction, both in MD and HRF-MDFT. In figure El we 
report the projection of P on the z axis (noted P z ) as a 
function of coordinate z {i.e. it is averaged in each plane 
z). Quantitative agreement is again found between the 
reference molecular dynamics simulations and the HRF- 
MDFT. A first peak is found at the location of the pre- 
peak, another one at the main peak, etc. Between each 
maximum, the sign of P z changes. Maps of P z in the 
planes of the prepcak, main peak and secondary peak 
are given in figure [7] Similarly to the density, the over- 
all shapes from both methods are in very good agree- 
ment. The polarization density is overestimated by HRF- 
MDFT in the prepeak, which is expected from previous 
work on small polarized systems^ where the polariza- 
tion close to the solute is often slightly overestimated. In 
the main and secondary peaks, the polarization densities 
arc found larger in MD than in HRF-MDFT, although 
not significantly. This may be due to a balance of the 
layer-by-layer polarizations in HRF-MDFT. 

In order to get more information on the orientational 
properties per molecule, we also define several orientation 
order parameters. One is the local orientation of the 
polarization density, defined as the cosine of the angle 
between the polarization density P(r) and the normal z 
to the surface, 



cos Op(r) 



l|P(r)ll 



Another is the averaged molecular orientation 



cosier) = 



n (r) 



(17) 



(18) 



We plot the average molecular orientation cosd^(z) = 
(cos 0/i(r)) in figure [3 The solvent molecules show 
strong average orientation in the first layer. On the con- 
trary, the high polarization density in the second, more 
cohesive layer reflects the large number of molecules in- 
side, each one with a relatively small preferential orien- 
tation along z. 

One can get more insight into the orientation proper- 
ties of the solvent molecules by looking at the order pa- 
rameter cos0p defined in Eq. [T7J When cosOp = 1 (—1), 
the dipolc is oriented against (toward) the surface on the 
left of the superccll, and vice versa for the other surface. 




20 30 
Z(A) 



40 



50 



FIG. 6. Component along the normal z to the surface of the 
dimensionless solvent polarization density, between two py- 
rophyllite layers, as calculated by molecular dynamics (open 
circles) and HRF-MDFT (black line). The density profile is 
also presented in arbitrary units (red dashed line). Maxima in 
the polarization density correspond to maxima in the solvent 
number density. 
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FIG. 7. Maps of the polarization density projected on the 
normal z to the surface axis, P z , in the planes of the pre- 
peak (left), main peak (center) and secondary peak (right), 
as calculated by molecular dynamics simulations (top) and 
HRF-MDFT (bottom). 



In figure [51 we plot cos 8p in the plane of the first three 
peaks identified earlier and in an intermediate coordinate 
z between the prepeak and the main peak. We observe 
that solvent molecules in the prepeak, strongly oriented, 
point against the surface. In the main peak, dipoles are 
preferentially oriented toward the surface when localized 
on top of O atoms. In the secondary layer, preferential 
polarization is found in the whole plane against the sur- 
face. 
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FIG. 8. Profile along the normal to the surface z of the av- 
erage molecular orientation, cos# M (z) (black line). The den- 
sity profile of the solvent is also given in arbitrary units (red 
dashed line). 



1 

0.5 


-0.5 
-1 









• 






1 A 








~ i 1 1 1 r 

b) 



6 9 12 15 18 

z(A) 



0.5 1 

Scale factor 



0.06 



0.04 



0.02 



FIG. 10. a) Number density profile along the z axis with 
point charges of the clay atoms scaled by a factor of 1 i.e. 
completely turned on (full black line), 0.5 (dashed red line), 
and i.e. turned off (blue dashed line). Only the prepeak 
is affected, b) Relative change in solvation free energy as a 
function of the scale factor. The relative changes stay below 
5.5 %. 



less than 6 % when charges are turned off. The contri- 
bution of the electrostatics is thus small compared to the 
Lennard- Jones one, which is not an obvious result even 
for a surface with a zero net charge. 



FIG. 9. From left to right, order parameter cos6>p in the 
prepeak, at the intermediate state between the prepeak and 
the main peak, in the main peak, and in the secondary peak. 



C. The Role of Electrostatics 

The computational efficiency of the molecular density 
functional theory unlocks systematic studies of the sol- 
vent structure and thermodynamics, such as the relative 
roles of van der Waals and electrostatic interactions to 
the CLAYFF force field. They are modeled respectively 
by Lennard- Jones interactions and a point charge distri- 
bution reported in table |TTJ Figure [TU] reports the num- 
ber density profile along the z axis for modified systems 
where the CLAYFF charges of the clay atoms have been 
scaled by a factor of 1, 0.5 and 0. It is observed that, sur- 
prisingly, only the shape of the prepeak is modified when 
the point charges are turned off. This part of the density 
profile evolves from a localized peak (with electrostatics 
on) to a shoulder of the main peak when quenched. The 
rest of the number density profile is unchanged. As might 
be expected, the polarization vanishes when charges are 
turned off. In figure [10] is also plotted the relative change 
in solvation free energy F [p] for scale factors between 
(turned off) and 1 (completely turned on) of the charges 
of the clay atoms. The solvation energy is affected by 



IV. CONCLUSION 

We have shown that the molecular density functional 
theory in the homogeneous reference fluid approximation 
(HRF-MDFT) is able to handle the study of solvation 
properties of a complex clay surface of several hundred 
atoms. The description of its structural and oricntational 
properties is quantitatively comparable to reference all- 
atom molecular dynamics, while reducing the CPU-time 
by two orders of magnitude. This means that the com- 
putation of these properties is now accessible on a simple 
workstations within minutes. The HRF-MDFT calcula- 
tions allowed to accurately describe the subtle structure 
of the first molecular layers at the surface, in particular 
the oricntational properties. The only noticeable differ- 
ence lies in a localized zone of high density at the center 
of hexagons, a defect due to a small overestimation of the 
polarization close to the solute inherent to the method. 
Some improvements in this respect are in progress. 

The numerical efficiency of our approach made it possi- 
ble to analyze the relative contributions of the CLAYFF 
force field to the solvation properties. The density pro- 
file and the solvation energy are found to be insensitive 
to the electrostatic contribution. This may imply that 
the role of electrostatics in charged clays may be rea- 
sonably reduced to their charged defects. We think that 
this finding may be of high interest considering the large 
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amount of work dedicated to the derivation of accurate 
point charge distributions when building force fields for 
such systems. 

Finally, the next step of this work will be to consider 
the solvation of clay surfaces by a realistic model of wa- 
ter, at either a dipolar or multipolar level, and by ionic 
solutions. 
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